//-L ../components_lib  -Wl,--gc-sections
#include <iostream>
#include <cctype>
#include <algorithm>
#include <vector>
#include <cassert>
#include <string>
#include <deque>
#include <iomanip>

// boost
#include <boost/tuple/tuple.hpp>
#include <boost/foreach.hpp>
#include <boost/regex.hpp>

// components
#include <print_sequence.h>
#include <print_error.h>
#include <cdna_prot.h>
using namespace cdna_prot;
#include <progress_indicator.h>
#include "get_arg.h"
#define _MY_STL_MACROS
#include <stlhelper.h>
#include <bcpp_parsefasta.h>
#include <md4.h>
#include <sequence_gaps.h>
#define DEBUG_LEVEL 9
#include <debug_func.h>
#include <commify.h>

using std::vector;
using std::deque;
using std::string;
using std::cerr;
using std::cout;
using std::sort;
using std::equal_range;
using std::istream;
using std::ostream;
using std::replace;
using std::remove;
using std::transform;
using std::reverse;

using boost::tie;
using boost::tuple;
using boost::make_tuple;
using BioCPP::get_seq_at_filepos;
using BioCPP::Process_FA_Sequences;
using BioCPP::Process_FA_Sequences_tellg;
//typedef regex::match_results::backref_type regex_bref;

namespace{
	inline unsigned to_uint(char c){return static_cast<unsigned>(c);}
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	index proteins

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
struct t_protid_filepos
{
	t_protid_filepos(const string& protid_, unsigned filepos_):
			protid(protid_), filepos(filepos_){}
	string		protid;
	unsigned	filepos;
	// order by id only
	bool operator < (const t_protid_filepos& o) const
	{
		return (protid.length() == o.protid.length() ?
				protid < o.protid	 :
				protid.length() < o.protid.length());

//		return protid <  o.protid;
	}
};

struct t_intron
{
	unsigned beg;
	unsigned len;
    t_intron(unsigned b, unsigned l):beg(b), len(l){	}
};

struct t_mismatch_aa
{
	unsigned beg;
	unsigned len;
	int intron;
	string aas;
	string aa_triframe;
    t_mismatch_aa(unsigned b, unsigned l, int intron_, const string& aas_, const string& aa_triframe_):
		beg(b), len(l), intron(intron_), aas(aas_), aa_triframe(aa_triframe_){}
};

typedef deque<t_intron> t_introns;
typedef deque<t_mismatch_aa> t_mismatch_aas;


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		find_consecutive_aa_matches

//			match peptide amino acid to cdna / virtual transcript for at least n codons
//			if successful, return true
//			cdna_pos and prot_pos will contain the positions at which the peptide/cdna
//				start to match

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
bool find_consecutive_aa_matches(	const string& tri_frame_aa, const string& prot_seq,
									unsigned& cdna_pos, unsigned& prot_pos,
									const unsigned consecutive_matching_aa)
{
	debug_func_cerr(5);
	// make sure there is enough cdna for a consecutive nth amino acid match
	unsigned cdna_limit	= tri_frame_aa.length() - (consecutive_matching_aa - 1) * 3;
	if (cdna_pos >= cdna_limit)
	{
		debug_cerr(5, "Not enough cDNA left");
		return false;
	}

	// savematch starting position for virtual transcript/cdna and peptide
	unsigned cdna_starting_pos = cdna_pos;
	unsigned prot_starting_pos = prot_pos;

	debug_cerr(6, "prot_pos=" << prot_pos);
	debug_cerr(6, "prot_seq.length()=" << prot_seq.length());
	debug_cerr(6, "cdna_pos=" << cdna_pos);
	debug_cerr(6, "cdna_limit=" << cdna_limit);
	debug_cerr_start(6, "Starting Match:");
	
	// make sure that we haven't run out of matching amino acids
	while (prot_pos + consecutive_matching_aa <= prot_seq.length())
	{
		// always start from this position along the virtual transcript / cdna
		cdna_pos   = cdna_starting_pos;
		while (cdna_pos < cdna_limit)
		{
			// be optimistic and see if nth next aa match
			unsigned unmatched_aa = consecutive_matching_aa;
			unsigned matching_pos = prot_pos;
			unsigned cdna_matching_pos = cdna_pos;
			unsigned prot_end = prot_seq.length();
			while (unmatched_aa && matching_pos != prot_end)
			{
				// do not match on gaps
				if (prot_seq[matching_pos] != '-')
				{
					// match
					if (tri_frame_aa[cdna_matching_pos] == prot_seq[matching_pos])
					{
						debug_run_block(6, if(unmatched_aa == consecutive_matching_aa)
											cerr << ", ";)
						debug_cerr_append(6, '+');
						cdna_matching_pos +=3;
						--unmatched_aa;
					}
					//
					//	some programmes, especially exonerate, put an X over a stop
					//
					else if (tri_frame_aa[cdna_matching_pos] == '*' && 
							 prot_seq[matching_pos] == 'X')
					{
						debug_run_block(6, if(unmatched_aa == consecutive_matching_aa)
											cerr << ", ";)
						debug_cerr_append(6, '*');
						cdna_matching_pos +=3;
						--unmatched_aa;
					}

					// no match
					else
						break;
				}
				++matching_pos;
			}

			// match
			if (!unmatched_aa)
			{
				debug_cerr(5, "==> Found contiguous match");
				return true;
			}

			// does not match, try shifting frames
			++cdna_pos;
		}

		// does not match at all with these amino acids
		// discard one and start from the beginning again
		++prot_pos;

	}

	// run out of amino acids to match
	// be good citizens and return variables to original values
	cdna_pos = cdna_starting_pos;
	prot_pos = prot_starting_pos;
	debug_cerr_append(6, '\n');
	debug_cerr(5, "Out of amino acids");
	return false;
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		map_cdna_to_prot

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
bool save_mismatch_cdna(t_introns& 		introns,
						unsigned		mismatch_cdna_pos,
						unsigned		cdna_pos)
{
	debug_func_cerr(5);
	if (mismatch_cdna_pos < cdna_pos)
	{
		// save missing intron
		introns.push_back(t_intron(	mismatch_cdna_pos, 
									cdna_pos - mismatch_cdna_pos));
		return true;
	}
	return false;
}

unsigned save_mismatch_prot(t_mismatch_aas& 	unmatched_aa,
							const string&		prot_seq,
							unsigned			mismatch_prot_pos,
							unsigned			prot_pos,
							int					intron_index,
							t_introns& 			introns,
							const string&		tri_frame_aa)
{
	debug_func_cerr(5);
	assert(mismatch_prot_pos <= prot_pos);

	// nothing missing
	if (mismatch_prot_pos >= prot_pos)
	{
		debug_cerr(5, "No missing aa");
		return 0;
	}

	// only missing hyphens: ignore
	assert(prot_pos <= prot_seq.length());
	string prot = prot_seq.substr(mismatch_prot_pos, prot_pos - mismatch_prot_pos);
	prot.erase(std::remove(prot.begin(),  prot.end(),  '-'),  prot.end());
	if (!prot.size())
	{
		debug_cerr(5, "No missing aa");
		return 0;
	}

	if (intron_index >= 0)
	{
		assert(intron_index <= (int)introns.size());
		string tri_frame;
		if (introns[intron_index].beg < tri_frame_aa.length() && 
			(introns[intron_index].beg + introns[intron_index].len - 2) <= tri_frame_aa.length())
			tri_frame = tri_frame_aa.substr(introns[intron_index].beg, 
																	introns[intron_index].len - 2);
		unmatched_aa.push_back(t_mismatch_aa(	mismatch_prot_pos, 
												prot_pos - mismatch_prot_pos,
												intron_index,
												prot, tri_frame));
	}															
	else
		unmatched_aa.push_back(t_mismatch_aa(	mismatch_prot_pos, 
												prot_pos - mismatch_prot_pos,
												intron_index,
												prot,
												string()));
	debug_cerr(5, "Normal exit");
	return prot.size();
}

unsigned map_cdna_to_prot(	const string&	cdna_index,
							const string&	prot_seq,
							string& 		cdna_mapped ,
							string& 		prot_mapped,
							unsigned 		consecutive_matching_aa,
							t_introns& 		introns,
							t_mismatch_aas& unmatched_aa)
							//,char use_X_for_aa_mismatches)
{
	debug_func_cerr(5);
	introns.clear();

	// make three frame translation
	string tri_frame_aa;
	const char* DNA = cdna_index.c_str();
	for (unsigned i = 0; i < cdna_index.length()-2; ++i)
		tri_frame_aa += t_cdna_index_to_prot()(DNA + i);
	consecutive_matching_aa = std::min(consecutive_matching_aa, tri_frame_aa.length() / 3);

/*	std::cerr<< tri_frame_aa << "\n";
	string cdna_index_txt = cdna_index;
	std::transform(beg_to_end(cdna_index_txt), cdna_index_txt.begin(), DNA_num_to_text());
	cerr << cdna_index_txt << "\n";
	std::cerr<< prot_seq << "\n";
*/

	unsigned prot_pos = 0;
	unsigned cdna_pos = 0;
	// the number of amino acids which we failed to match

	unsigned cnt_unmatched_aa	= 0;
	unsigned mismatch_cdna_pos	= 0;
	unsigned mismatch_prot_pos	= 0;
	while (cdna_pos < tri_frame_aa.length()	&&
		   prot_pos < prot_seq.length())
	{
		//
		// insert aligning hyphens
		//
		if (prot_seq[prot_pos] == '-')
		{
			prot_mapped += '-';
			cdna_mapped.append(3, '\x00');
			++prot_pos;
			continue;
		}

		//
		// matched: increment and continue
		//
		if (tri_frame_aa[cdna_pos] == prot_seq[prot_pos])
		{
			prot_mapped += prot_seq[prot_pos];
			cdna_mapped += cdna_index.substr(cdna_pos, 3);
			cdna_pos += 3;
			prot_pos++;
			mismatch_cdna_pos	= cdna_pos;
			mismatch_prot_pos	= prot_pos;
			continue;
		}

		//
		//	some programmes, especially exonerate, put an X over a stop
		//
		if (tri_frame_aa[cdna_pos] == '*' && prot_seq[prot_pos] == 'X')
		{
			prot_mapped += 'X';
			cdna_mapped.append(3, '\x05');	// N
			cdna_pos += 3;
			prot_pos++;
			mismatch_cdna_pos	= cdna_pos;
			mismatch_prot_pos	= prot_pos;
			continue;
		}

		//
		// Here, must be mismatch. Save position
		mismatch_cdna_pos	= cdna_pos;
		mismatch_prot_pos	= prot_pos;


		//
		//	Initial match must be > threshold length
		//  	failed: mismatch from start
		//
		if (prot_mapped.length() < consecutive_matching_aa)
		{
			debug_cerr(5, "Beginning too short");
			prot_mapped.clear();
			cdna_mapped.clear();
			mismatch_cdna_pos   = 
			mismatch_prot_pos   = 
			prot_pos			= 
			cdna_pos			= 0;
		}



		//
		//	Try to find next consecutive run of matching cDNA/aa
		//
		if (! find_consecutive_aa_matches(tri_frame_aa, prot_seq, cdna_pos, prot_pos,
											consecutive_matching_aa))
			break;
		debug_cerr(5, "Resumed at " << prot_pos << "," << cdna_pos
				    << "[" << prot_mapped << "]");

		// special handling for cases where one or more 'N's are missing
		char	 mismatch_aa		= prot_seq[mismatch_prot_pos];
		if (mismatch_aa == 'X' && mismatch_prot_pos < prot_pos)
		{
			if (cdna_pos > mismatch_cdna_pos)
			{
				prot_mapped += 'X';
				mismatch_prot_pos++;

				if (cdna_pos >= mismatch_cdna_pos + 2)
					cdna_mapped += '\x05' + cdna_index.substr(mismatch_cdna_pos, 2);
				else if (cdna_pos == mismatch_cdna_pos + 1)
					cdna_mapped += string(2, '\x05') + cdna_index.substr(mismatch_cdna_pos, 1);
				else
				{
					assert (cdna_pos == mismatch_cdna_pos);
					cdna_mapped.append(3, '\x05');
				}
			}
		}

		// codon may be split across splice site
		else if (mismatch_prot_pos && 
				 mismatch_prot_pos < prot_pos && cdna_pos - mismatch_cdna_pos >= 4)
		{
			string spliced_codon = 	cdna_index[mismatch_cdna_pos] +
									cdna_index.substr(cdna_pos - 2, 2);
			if (t_cdna_index_to_prot()(spliced_codon.c_str()) == mismatch_aa)
			{
				cdna_mapped += spliced_codon;
				prot_mapped += mismatch_aa;
				mismatch_prot_pos++;
			}
			else
			{
				spliced_codon = 	cdna_index.substr(mismatch_cdna_pos, 2) +
									cdna_index[cdna_pos - 1];
				if (t_cdna_index_to_prot()(spliced_codon.c_str()) == mismatch_aa)
				{
					cdna_mapped += spliced_codon;
					prot_mapped += mismatch_aa;
					mismatch_prot_pos++;
				}
			}
		}


		// add the requisite number of 'X's if necessary to fill in missing residues
		if (mismatch_prot_pos < prot_pos)
			// && use_X_for_aa_mismatches)
		{
#ifndef NDEBUG
			unsigned cnt_xx = 0;
#endif
			for (unsigned i = mismatch_prot_pos; i < prot_pos; ++i)
			{
				if (prot_seq[i] == '-')
				{
					prot_mapped += '-';
					cdna_mapped.append(3, '\x00');
				}
				else
				{
#ifndef NDEBUG
					++cnt_xx;
#endif
					//use_X_for_aa_mismatches);
					prot_mapped += '\x01';
					cdna_mapped.append(3, '\x06');
				}
			}
			debug_cerr(5, "added "	<< cnt_xx
									<< "(between " << mismatch_prot_pos << " and " << prot_pos
									<< ") Xs");
		}


		//
		//	mis-matching amino acid / cdna
		//
		debug_cerr(5, "save mismatches " << (cdna_pos - mismatch_cdna_pos) << ":"
										 << (prot_pos - mismatch_prot_pos));
		bool is_intron = save_mismatch_cdna(introns, mismatch_cdna_pos, cdna_pos);
		cnt_unmatched_aa += save_mismatch_prot(unmatched_aa, prot_seq, mismatch_prot_pos, prot_pos, 
											  is_intron ? int(introns.size() - 1): -1,
											  introns, tri_frame_aa);
		mismatch_cdna_pos	= cdna_pos;
		mismatch_prot_pos	= prot_pos;
	}

	// 
	// save 3' utr
	// 
	while (prot_pos < prot_seq.length())
	{
		if (prot_seq[prot_pos] == '-')
		{
			prot_mapped += '-';
			cdna_mapped.append(3, '\x00'); // hyphens
		}
		else // if (use_X_for_aa_mismatches)
		{
			cdna_mapped.append(3, '\x06');	// N
			prot_mapped += '\x01';			// use_X_for_aa_mismatches;
		}
		++prot_pos;
	}
	cdna_pos = cdna_index.length();

	bool is_intron = save_mismatch_cdna(introns, mismatch_cdna_pos, cdna_pos);
	cnt_unmatched_aa += save_mismatch_prot(unmatched_aa, prot_seq, mismatch_prot_pos, prot_pos, 
										  is_intron ? int(introns.size() - 1): -1,
											  introns, tri_frame_aa);

	debug_cerr(5, "Normal exit");

	// return the number of amino acids which we failed to match
	return cnt_unmatched_aa;

}


struct t_index_protid
{
	t_index_protid(	vector<t_protid_filepos>&		map_protid_filepos_,
					const string&					regex_str_,
					ostream&						verbose_strm_):

					map_protid_filepos (map_protid_filepos_),
					regex_str(regex_str_),
					count_invalid_seq(0),
					count_no_regex_match(0),
					progress(verbose_strm_, 1000),
					//don't allow hyphens here because any sequence will consists
					// only of hyphens is seriously doggy
					valid_seq_char("ACDEFGHIKLMNPQRSTVWXY", true),
					verbose_strm(verbose_strm_)
		 {}

	bool operator()(const string& accession, string sequence,
					unsigned acc_filepos, unsigned seq_filepos);

	vector<t_protid_filepos>& 	map_protid_filepos;
	string						regex_str;
	unsigned					count_invalid_seq;
	unsigned					count_no_regex_match;

	t_progress_indicator		progress;
	cdna_prot::t_validate_prot	valid_seq_char;
	ostream& 					VERBOSE(){return verbose_strm;}
	ostream&					verbose_strm;
};


bool t_index_protid::operator()(	const string& accession,
									string	sequence,
									unsigned acc_filepos,
									unsigned )
{

	// get rid of line breaks, extraneous characters etc.
	valid_seq_char.validate(sequence);
	if (!sequence.length())
		++count_invalid_seq;
	else
	{
		// if regex does not have bracket, get whole string
		static const boost::regex acc_code_regex(regex_str);
		static unsigned regex_use_bracket = (regex_str.find('(') == string::npos) ? 0 : 1;
		boost::match_results<string::const_iterator> acc_code_match;

		if (!regex_search(accession, acc_code_match, acc_code_regex) ||
			 acc_code_match[regex_use_bracket].length() == 0)
			++count_no_regex_match;
		else
		{
			map_protid_filepos.push_back(
						t_protid_filepos(	acc_code_match[regex_use_bracket],
											acc_filepos));
		}

	}

	++progress;
	return true;
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	process cdna

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
struct t_align_cdna_to_prot
{
	t_align_cdna_to_prot(vector<t_protid_filepos>&	map_protid_filepos_,
						 const string&				regex_str_,
						 unsigned					consecutive_matching_aa_,
						 istream&					f_prot__in_,
						 ostream&					f_prot_out_,
						 ostream&					f_cdna_out_,
						 ostream&					f_intron_out_,
						 ostream&					f_cdna_tout_,
						 ostream&					f_missing_identifiers_,
						 bool						use_f_err__log_,
						 ostream&					f_err__log_,
						 char						use_X_for_aa_mismatches_,
						 ostream&					verbose_strm_):

		map_protid_filepos (map_protid_filepos_),
		regex_str(regex_str_),
		count_invalid_seq(0),
		count_no_regex_match(0),
		count_no_matching_prot(0),
		count_mapped(0),
		count_aa_mismatches(0),
		progress(verbose_strm_, 1000),
		f_prot__in(f_prot__in_),
		f_prot_out(f_prot_out_),
		f_cdna_out(f_cdna_out_),
		f_intron_out(f_intron_out_),
		f_cdna_tout(f_cdna_tout_),
		f_missing_identifiers(f_missing_identifiers_),
		use_f_err__log(use_f_err__log_),
		f_err__log(f_err__log_),
		consecutive_matching_aa(consecutive_matching_aa_),
		use_X_for_aa_mismatches(use_X_for_aa_mismatches_),
		verbose_strm(verbose_strm_),
		// allow for aligned sequences
		valid_seq_char("ACDEFGHIKLMNPQRSTVWXY-", true)
		 {}


	bool operator()(const string& accession, string sequence,
					unsigned line_cnt);
	ostream& 					VERBOSE(){return verbose_strm;}

	vector<t_protid_filepos>& 	map_protid_filepos;
	string						regex_str;
	unsigned					count_invalid_seq;
	unsigned					count_no_regex_match;
	unsigned					count_no_matching_prot;
	unsigned					count_mapped;
	unsigned					count_aa_mismatches;

	t_progress_indicator		progress;
	istream&					f_prot__in;
	ostream&					f_prot_out;
	ostream&					f_cdna_out;
	ostream&					f_intron_out;
	ostream&					f_cdna_tout;
	ostream&					f_missing_identifiers;
	bool						use_f_err__log;
	ostream&					f_err__log;
	unsigned					consecutive_matching_aa;
	char						use_X_for_aa_mismatches;
	ostream&					verbose_strm;
	cdna_prot::t_validate_prot	valid_seq_char;
};


void print_tri_frame(ostream& os, const string& str)
{
	for (unsigned offset = 0; offset < 3; ++offset)
	{
		os << "frame " << (offset + 1) << ": ";
		os << string(offset + 1, ' ');
		for (unsigned i = offset; i < str.length(); i += 3)
			os << str[i] << "  ";
		os << "\n";
	}
}
bool t_align_cdna_to_prot::operator()(	const string& cdna_acc, string cdna_index,
										unsigned )
{
	debug_func_cerr(5);
	debug_cerr(5, cdna_acc);
	++progress;




	//
	// get cdna protein id
	//
	// if regex does not have bracket, get whole string
	static const boost::regex acc_code_regex(regex_str);
	static unsigned regex_use_bracket = (regex_str.find('(') == string::npos) ? 0 : 1;
	boost::match_results<string::const_iterator> acc_code_match;

	if (!regex_search(cdna_acc, acc_code_match, acc_code_regex) ||
		 acc_code_match[regex_use_bracket].length() == 0)
	{
		++count_no_regex_match;
		return true;
	}





	//
	// get matching protein sequence filepos
	//
	t_protid_filepos protid (acc_code_match[regex_use_bracket], 0);
	vector<t_protid_filepos>::iterator i, j;
	boost::tie(i,j) =	equal_range(beg_to_end(map_protid_filepos), protid);
	if (i == j)
	{
		f_missing_identifiers
						<< string(acc_code_match[regex_use_bracket])
						<< "\t(from the cdna file) was not found in the list of protein ids.\n";
		++count_no_matching_prot;
		return true;
	}
	unsigned prot_filepos = i->filepos;



	//
	// get matching protein sequence
	//
	string prot_acc;
	string prot_seq;
	f_prot__in.clear();
	get_seq_at_filepos(f_prot__in, prot_filepos, prot_acc, prot_seq);
	string orig_seq = prot_seq;
	valid_seq_char.validate(prot_seq);



	//
	// make sure cdna sequence is valid and convert to 1..6 index
	//
	// 		pad with extra Ns just in case because of degenerate genetic code.
	cdna_index += 'N';
	validate_cdna_to_index(cdna_index);
	if (!cdna_index.length())
	{
		++count_invalid_seq;
		return true;
	}

	// align cdna and proteins sequence
	string cdna_index_mapped, prot_seq_mapped;
	t_introns 		introns;
	t_mismatch_aas	unmatched_aa;

	unsigned cnt_aa_mismatches =	0;
	// try complete match first
	if (map_cdna_to_prot(	cdna_index, 	prot_seq,
							cdna_index_mapped, 	prot_seq_mapped,
							prot_seq.length() - 1,
							introns,
							unmatched_aa))
	{
		introns.clear();
		unmatched_aa.clear();
		cdna_index_mapped =
		prot_seq_mapped = string();
		cnt_aa_mismatches =	map_cdna_to_prot(	cdna_index, 	prot_seq,
															cdna_index_mapped, 	prot_seq_mapped,
															consecutive_matching_aa,
															introns,
															unmatched_aa);
	}
	
	debug_cerr(5, cnt_aa_mismatches);
	
	// if too many mistakes, try 3' to 5'
	if (cnt_aa_mismatches > 5)
	{
		string cdna_index_mapped_3prime, prot_seq_mapped_3prime;
		t_introns 		introns_3prime;
		t_mismatch_aas	unmatched_aa_3prime;
		reverse(beg_to_end(cdna_index));
		complement_cdna_index(cdna_index);
		// pad with extra Ns just in case because of degenerate genetic code.
		cdna_index += 0x5;
		debug_cerr(5, "Try matching in the 3' direction:");
		unsigned cnt_aa_mismatches_3prime =	map_cdna_to_prot(cdna_index, prot_seq,
														 	 cdna_index_mapped_3prime,
															 prot_seq_mapped_3prime,
															 consecutive_matching_aa,
															 introns_3prime,
															 unmatched_aa_3prime
															 );
															 //use_X_for_aa_mismatches);
		// if 3' more successful, use that instead
		if (cnt_aa_mismatches > cnt_aa_mismatches_3prime)
		{
			cnt_aa_mismatches	= cnt_aa_mismatches_3prime;
			cdna_index_mapped	= cdna_index_mapped_3prime;
			prot_seq_mapped		= prot_seq_mapped_3prime;
			introns.swap(introns_3prime);
			unmatched_aa.swap(unmatched_aa_3prime);
		}
		else
		{
			cdna_index.erase(cdna_index.length() - 1);
			reverse(beg_to_end(cdna_index));
			complement_cdna_index(cdna_index);
		}
	}
	assert(cdna_index_mapped.length() == prot_seq_mapped.length() *3);


	string cdna_seq_mapped = cdna_index_to_text(cdna_index_mapped);

	// record errors
	if (cnt_aa_mismatches)
	{
		debug_cerr(5, "Record errors");
		++count_aa_mismatches;
		if (use_f_err__log)
		{
			f_err__log
						<<  "\n!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
						<< ">" << protid.protid << "\n"
						<< "length=\t" << prot_seq.length()
						<< "\tmismatches=\t" << cnt_aa_mismatches << "\n";
			orig_seq.erase(remove(beg_to_end(orig_seq), '\n'), orig_seq.end());

			// show missing as '^'
			string prot_seq_mapped_tmp = prot_seq_mapped;
			replace(prot_seq_mapped_tmp.begin(), prot_seq_mapped_tmp.end(), '\x01', '^');
			f_err__log
				<< "original  " << orig_seq << "\n"
				<< "validated " << prot_seq << "\n"
				<< "mapped    " << prot_seq_mapped_tmp << "\n";
		}

		BOOST_FOREACH(t_mismatch_aa& wrong_aa, unmatched_aa) 
		{
			f_err__log	
						<<  "-------------------------------\n"
							"[Missing residues from position\t"
						<< (wrong_aa.beg + 1)
						<< "\tlength=\t" << wrong_aa.len << "]\n"
						"Missing:  ";

			BOOST_FOREACH(char c, wrong_aa.aas) 
				f_err__log << c << "  ";
			f_err__log << "\n";
			if (wrong_aa.intron >= 0)
			{
				string orig_cdna = cdna_index_to_text(cdna_index);
				print_tri_frame(f_err__log, wrong_aa.aa_triframe);
				f_err__log	<< "cDNA   : "
							<< orig_cdna.substr(introns[wrong_aa.intron].beg, 
												introns[wrong_aa.intron].len) << "\n";
				f_err__log << "Intron#: " << wrong_aa.intron << "\n";
			}
		}
	}

	{
		debug_cerr(5, "Introns");
		f_intron_out << ">" << protid.protid << "\n";
		BOOST_FOREACH(t_intron& intron, introns) 
		{
			debug_cerr(5, "!"  << intron.beg << "\t" << intron.len << " out of " << cdna_index.length());
			f_intron_out << intron.beg << "\t" << intron.len << "\n";
#ifndef	NDEBUG
			// make three frame translation
			string tri_frame_aa;
			string orig_cdna = cdna_index_to_text(cdna_index);
			const char* DNA = cdna_index.c_str() + intron.beg;
			// intron len - 2 could be less than 2 if intron is less than a codon
			if (intron.len >= 3)
			{
				for (unsigned i = 0; i < intron.len - 2; ++i)
				{
					tri_frame_aa += t_cdna_index_to_prot()(DNA + i);
				}
				print_tri_frame(f_intron_out, tri_frame_aa);
				f_intron_out << orig_cdna.substr(intron.beg,  intron.len) << "\n";
			}

#endif
		}
		debug_cerr(5, "Finished Introns");

	}


	//
	// don't do anything if failed to match
	//
	if (!prot_seq_mapped.length())
		return true;

	//
	// replace or delete missing residues
	//
	debug_cerr(5, "Replace or delete missing residues");
	if (use_X_for_aa_mismatches)
	{
		replace(prot_seq_mapped.begin(), prot_seq_mapped.end(), '\x01', use_X_for_aa_mismatches);
		replace(cdna_seq_mapped.begin(), cdna_seq_mapped.end(), '^', 'N');
	}
	else
	{
		prot_seq_mapped.erase(remove(prot_seq_mapped.begin(), prot_seq_mapped.end(), '\x01'), prot_seq_mapped.end());
		cdna_seq_mapped.erase(remove(cdna_seq_mapped.begin(), cdna_seq_mapped.end(), '^'), cdna_seq_mapped.end());
	}

	//
	// nothing to align: fail
	//
	if (!cdna_seq_mapped.length())
		return true;

	t_digest seq_digest;


	// output cdna as fasta
	debug_cerr(5, "Output");
	print_fasta(f_cdna_out, cdna_acc, cdna_seq_mapped);

	// output pi, prot_id, cdna sequence as a table
	f_cdna_tout		<< md4_hash(prot_seq_mapped, seq_digest).str('\t')  << '\t'
					<< protid.protid									<< '\t'
					<< cdna_seq_mapped									<< '\t'
					<< md4_hash(cdna_seq_mapped, seq_digest).str('_')   << '\n';


	// output protein sequence in fasta format
	print_fasta(f_prot_out, prot_acc, prot_seq_mapped);


	count_mapped++;
	return true;
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
/*
43738   protein sequences parsed.
43738   cdna sequences processed.
43738   cdna sequences mapped successfully.
16      peptide sequences having 1 or more residues without matching codons.
*/

void summary(t_program_args& args, int count, const string& message, bool log_error = false)
{
	args.VERBOSE()	<< "\t"
					<< std::setw(7) << commify(count)
					<< " " + message + "\n";
	args.SUMMARY()
					<< std::setw(7) << commify(count)
					<< " " + message + "\n";
	if (log_error)
	{
		args.CERR()	<< "\t"
						<< std::setw(7) << commify(count)
						<< " " + message + "\n";
	}
}

int main (int argc, char *argv[])
{
	debug_func_cerr(5);
	// get command line arguments

	// if no arguments print usage and exit
	if (argc == 1)
	{
		char * my_argv[2] ={ argv[0], "-h"};
		int argc = 2;
		t_program_args args;
		if (!process_command_line_options(argc, my_argv, args))
			return 1;
	}
	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;




	//
	// index proteins
	//
	args.VERBOSE()	<<  "\n"
					<<	"\tIndexing list of protein sequences from [ "
					<<	args.prot__in	<<	" ] ...\n";

	vector<t_protid_filepos> map_protid_filepos;


	t_index_protid index_of_proteins(map_protid_filepos, args.regex_str, args.VERBOSE());
	Process_FA_Sequences_tellg(args.f_prot__in, index_of_proteins);
	index_of_proteins.progress.finish();
	sort(beg_to_end(map_protid_filepos));

	summary(args, map_protid_filepos.size(), "protein sequences parsed.");

	if (index_of_proteins.count_invalid_seq)
		summary(args, index_of_proteins.count_invalid_seq, "invalid protein sequences parsed.");

	if (index_of_proteins.count_no_regex_match)
		summary(args, index_of_proteins.count_no_regex_match, "protein sequences without protein IDs.");


	//
	// process cdna
	//
	args.VERBOSE()	<<	"\tAligning cdna from [ "
					<<	args.cdna__in	<<	" ] to matching protein sequences...\n";

	// output pi, prot_id, cdna sequence
	t_align_cdna_to_prot align_cdna_to_prot(	map_protid_filepos,
												args.regex_str,
												args.consecutive_matching_aa,
												args.f_prot__in,
												args.f_prot_out,
												args.f_cdna_out,
												args.f_intron_out,
												args.f_cdna_tout,
												args.f_missing_identifiers,
												args.logging_errors(),
												args.CERR(),
												args.use_X_for_aa_mismatches,
												args.VERBOSE());
	Process_FA_Sequences(args.f_cdna__in, align_cdna_to_prot);
	unsigned cnt_cdna = align_cdna_to_prot.progress.value();
	align_cdna_to_prot.progress.finish();


	//
	// statistics
	//
	summary(args, cnt_cdna, "cdna sequences processed.");
	summary(args, align_cdna_to_prot.count_mapped, "cdna sequences mapped successfully.");

	if (align_cdna_to_prot.count_invalid_seq)
		summary(args, align_cdna_to_prot.count_invalid_seq, 
					"invalid cdna sequences parsed.", 
					args.logging_errors());

	if (align_cdna_to_prot.count_no_regex_match)
		summary(args, align_cdna_to_prot.count_no_regex_match, 
				"cdna sequences without protein IDs.", 
				args.logging_errors());

	if (align_cdna_to_prot.count_no_matching_prot)
		summary(args, align_cdna_to_prot.count_no_matching_prot, 
				"cdna sequences without corresponding proteins IDs.", 
				args.logging_errors());

	if (align_cdna_to_prot.count_aa_mismatches)
		summary(args, align_cdna_to_prot.count_aa_mismatches, 
				"peptide sequences having 1 or more residues "
				"without matching codons.", 
				args.logging_errors());


	args.VERBOSE()	<< "\n";





/*

	print STDERR    "\t", "-" x 6, "\n";
    print STDERR	"\t", $cnt_success, "\tcdna mapped onto their corresponding proteins.\n";

	print "*"x80, "\n";



  */


	}
	catch(std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}


	return(0);
}


